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The expected phenomenology of non-interacting topological band insulators (TBI) is now largely 
theoretically understood. However, the fate of TBIs in the presence of interactions remains an 
active area of research with novel, interaction-driven topological states possible, as well as new 
exotic magnetic states. In this work we study the magnetic phases of an exchange Hamiltonian 
arising in the strong interaction limit of a Hubbard model on the honeycomb lattice whose non- 
interacting limit is a two-dimensional TBI recently proposed for the layered heavy transition metal 
oxide compound, (Li,Na)2lr03. By a combination of analytical methods and exact diagonalization 
studies on finite size clusters, we map out the magnetic phase diagram of the model. We find that 
strong spin-orbit coupling can lead to a phase transition from an antiferromagnetic Neel state to a 
spiral or stripy ordered state. We also discuss the conditions under which a quantum spin liquid may 
appear in our model, and we compare our results with the different but related Kitaev-Heisenberg- 
J2-J3 model which has recently been studied in a similar context. 

PACS numbers: 71.10.Fd,71.10.Pm,73.20.-r 



I. INTRODUCTION 

Topological band insulators (TBI) preserve time re- 
versal symmetry, have a bulk band gap originating in 
strong spin-orbit coupling, and possess the physical sig- 
nature of an odd number of gapless Dirac nodes at time- 
reversal invariant momenta on boundary. 1-7 The spin- 
momentum locking originating from strong spin-orbit 
coupling in these systems and odd number of Dirac points 
results in gapless boundary states immune to Ander- 
son localization. 8 Since the experimental discovery 9 ' 10 
of these topological states of matter, there has been an 
explosion of both theoretical 11-14 and experimental 15-20 
works aimed at unraveling the fascinating properties of 
topological insulators. Noninteracting topological insu- 
lators and superconductors can be classified in terms of 
topological invariants 2,4,21 ' 22 where the band topology 
fully characterizes the properties of the insulator and su- 
perconductor. What remains to be understood is the 
full effect of Coulomb interactions, including the possi- 
ble magnetic phases which could arise from the nontrivial 
band topology in the limit of intermediate to strong inter- 
actions. Almost all experimentally discovered topologi- 
cal insulators to date can be understood within a single- 
particle picture (i.e., band theory) where the correlation 
effects are weak. 23-26 Besides the known topological insu- 
lators, there are scores of other materials whose noninter- 
acting and weakly interacting limits have been predicted 
to be topological insulators. 11-13,27 

It is well known that the nontrivial band topology orig- 
inates from relativistic spin-orbit coupling strong enough 
to cause a band inversion at an odd number of time- 
reversal invariant momenta in the Brillouin zone. 3 ' 28 
Hence, it is natural to expect that the materials with 
heavy elements with their large spin-orbit coupling may 



provide fertile ground in the search for topological insu- 
lators. Among them, the transition metal oxides with 5d 
elements such as Ir and Os attract attention as the 5d 
orbitals are spatially more extended, and therefore less 
correlated compared to the 3d and 4d orbitals. A few ex- 
amples are the pyrochlore iridates, A2IX2O7 (A is a rare 
earth element), the layered compounds (Li, Na)2lr03 and 
Sr 2 Ir04, and the hyperkagome Na 2 Ir,i08. Naively, these 
materials are expected to be metallic with a partially 
filled transition metal ion shell and a relatively small ef- 
fective Hubbard U (on the order of 0.5-2.0 eV). However, 
since the intrinsic spin-orbit coupling is strong (on the 
order of 0.4-1.0 eV), the spin and orbital degrees of free- 
dom are entangled 29 which can dramatically influence the 
band topology of the systems. 30-33 

The effect of electron interactions on band topology 
has been studied for variety of models. A particularly 
well studied model is the famous Kane-Mele model 1 ' 2 
with a Hubbard interaction added. 34-41 In the absence 
of spin-orbit coupling, the model reduces to the Hub- 
bard model on the honeycomb lattice. Quantum Monte 
Carlo simulations were used to map out the phase dia- 
gram of the Hubbard model. Three phases were found: 
(i) a metallic phase, (ii) a quantum spin liquid (QSL), and 
(iii) an antiferromagnetic (AFM) phase with increasing 
the strength of local Hubbard term. 42 While the metal- 
lic phase is immediately converted into the quantum spin 
Hall state (QSH) upon the inclusion of a second-neighbor 
spin-orbit coupling, the QSL phase is stable over a small 
range of spin-orbit coupling. 38-41 Both QSH and QSL 
phases are unstable in the strong interacting limit where 
an in-plane antiferromagnetic ordering arises. 36 ' 38,43 

Besides the two-dimensional Kane-Mele-Hubbard 
model, three dimensional systems have drawn attention. 
Pyrochlore oxides with heavy transition elements have 
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been studied using a strong spin-orbit coupling approach 
that splits the ti g manifold into a lower j = 3/2 man- 
ifold and a higher j = 1/2 manifold. The latter acts 
effectively as a spin-1/2 degree of freedom. When slave- 
particle approaches are applied, exotic topological Mott 
insulators with topologically protected gapless bound- 
ary spin excitations appear in a range of intermedi- 
ate strength Hubbard interactions. 31,32,44 If time-reversal 
symmetry is broken, Weyl semi-metals may also appear 
in the pyrochlores 45-48 and possibly also an axion insu- 
lator phase proximate to the Weyl semi-metal. 45,47,49 

In addition to the studies mentioned above, there 
are other works addressing the physics of interaction- 
generated spin-orbit coupling which could drive the sys- 
tem to a phase with nontrivial band topology. In this 
case, the topological order appears via spontaneously 
generated complex hopping terms which mimic those of 
an intrinsic spin-orbit coupling. For two dimensional sys- 
tems with quadratic band touching points in their non- 
interacting band structure, the leading instability would 
be a quantum anomalous Hall (QAH) effect and/or a 
topological insulator which, respectively, break the time 
reversal symmetry and spin rotational symmetry. 50-55 
Similar physics is also believed to possibly generate topo- 
logical phases in transition metal oxide heterostructures 
derived from the much lighter 3d elements. 48,56,57 

In the strong coupling limit, spin-orbit coupling can 
also affect the magnetic phases of the transition metal 
oxides. In the absence of spin-orbit coupling and orbital 
degeneracy the strong coupling limit can often be ade- 
quately described in terms of a pure spin Hamiltonian 
of the Heisenberg form. This is believed to be the case 
in the insulating parent compound of cuprate supercon- 
ductors, for example. However, an orbital degeneracy is 
often present in transition metal ions leading to a spin- 
orbital exchange interaction. 58 In contrast to a single 
orbital model, the resulting exchange interaction could 
be highly anisotropic and frustrated. 59,60 The entangled 
spin and orbital states break the SU(2) symmetry of the 
magnetic Hamiltonian giving rise to realizations of exotic 
spin models such as the Kitaev 61 or Heisenberg-Kitaev 
models 62,63 in transition metal oxides. 

Recently, the layered perovskites (Li, Na)2lr03 have 
been suggested to host exotic phases. Temperature de- 
pendent electrical resistivity and magnetic measurements 
clearly indicate their insulating nature with enhanced 
magnetic correlations at low temperatures. 64 The insu- 
lating ground state is thought to be interaction-driven 
and magnetically ordered at low temperatures. 65 Al- 
though strong Coulomb interactions make the realization 
of topological band insulators unlikely, the intrinsic spin- 
orbit coupling does substantially modify the effective spin 
Hamiltonian: The Heisenberg-Kitaev model 63 has been 
proposed to explain the strong suppression of magnetic 
correlations due to the possible proximity to a quantum 
spin liquid phase (T/v rj 15K is much smaller than the 
Curie-Weiss temperature 6 = — 116 66 ), though recent 
x-ray magnetic scattering experiments suggest the sys- 



tem is magnetically zig-zag ordered. 6 Subsequent works 
based on magnetic models might be able to describe this 
ordered phase. 68-72 

In this work, we consider an alternative magnetic 
Hamiltonian that is obtained in the strong interaction 
limit of the model introduced in Ref. [30], which at the 
nonintcracting level exhibits a two-dimensional Z 2 topo- 
logical insulator, and in the intermediate interaction 
regime may exhibit a novel interaction-driven topological 
insulator with a non-trivial ground state degeneracy and 
topologically protected collective modes. 73 In our deriva- 
tion of the effective spin Hamiltonian, we explicitly in- 
clude the second-neighbor real hopping (in addition to 
the complex hopping) predicted from band theory which 
then gives rise to an anisotropic exchange coupling, which 
fully breaks the SU(2) spin symmetry [see Eq.(3)]. By a 
combination of analytical and exact diagonalization stud- 
ies we map out the phase diagram of the model. The full 
phase diagram is shown in Fig. 8 where a variety of mag- 
netic phases result from the interplay between spin-orbit 
coupling and correlations. We also discuss the partial rel- 
evance of the model to magnetic ordered state discovered 
in layered (Li, Na) 2 Ir03 as the ground state with stripy 
order still possess some degree of zig-zag ordering. 

Our paper is organized as follows. In Sec. II we in- 
troduce both the non-interacting and magnetic exchange 
Hamiltonians. In Sec. Ill we use Schwinger Boson mean 
field theory (SBMFT) to address the ineffectiveness of 
the anisotropic exchange term in stabilizing a spin liquid 
phase. We then study possible classical magnetic phases 
in Sec. IV, and in Sec.V exact diagonalization is used to 
study the various magnetic phases and phase transitions 
between them. We further study the locations of criti- 
cal points by use of fidelity in Sec. VI, and present our 
conclusions in Sec. VII. Some details of the SBMFT are 
included in Appendix A and B. 



II. MODEL AND METHODS 

The layered oxides (Li, Na)2lr03 are composed of 
NaIrC"3 layers stacked along the c-axis and separated by 
a layer of Na [see Fig. 1(a)] (and likewise for the Li-based 
material), where the transition metal ions Ir +4 (5<i 5 ) site 
on the vertices of a honeycomb lattice. The measured 
magnetic moment is /i e // = 1-8 f-B verifying the de- 
scription of the model in terms of local moments with 
S e ff = 1/2. 64 The noninteracting model was argued, 
from a tight-binding fit to a density functional theory cal- 
culation, to be described by the following Hamiltonian, 30 

H = -t ]T 4 a c ia + I>l«*a^ (!) 

<ij>a -Cij3>7 afi 

where c\(ca) stands for the creation (annihilation) oper- 
ator of an electron in a spin-orbital coupled pseudospin- 
1/2. The first term describes spin- independent nearest- 
neighbor hopping, while the second term describes 
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(a) (b) 

FIG. 1. (Color online) (a) Lattice representation of Na2lr03. 
The blue (grey), red (dark) and yellow (light) balls indicate 
Iridium, Oxygen and Sodium, respectively. The Ir +4 ions are 
located on the vertices of the honeycomb lattice of NalrOg, 
stacked along c-axis. 74 (b) A top view of a honeycomb lat- 
tice. Solid lines stand nearest-neighbor (NN) coupling and red 
(dotted-dashed), green (dotted) and blue (dashed links) indi- 
cate the spin-dependent next-nearest-neighbor (NNN) hop- 
ping terms in Eq.(l) involving the Pauli spin matricies a x , 
a v , and a z , respectively. The same NNN lines also stand for 
the anisotropic exchange couplings in Eq.(3). Vertices shown 
by light and blue circles, squares and triangles indicate the 
local transformations on the triangular lattice which restores 
the SU(2) symmetry of third term in Eq.(3) [See text around 
Eq.(ll)]. 



the limit of strong Hubbard interaction (at half-filling in 
the j = 1/2 band) which results in the exchange Hamil- 
tonian 

// Y^.I^S.-S, .1, ]T -S,-- 2S7S]), (3) 

where J l3 =Ji = At 2 /U and J tj =J 2 = A{f) 2 /U stands for 
first and second neighbor links, respectively, [both con- 
tained in the first term of (3)], and J 3 — A(t") 2 /U denotes 
the strength second neighbor contributions coming from 
the imaginary spin-dependent hopping term. The Si de- 
notes the effective spin- 1/2 moment of the Ir 4+ ions, and 
the 7 indicates the direction of the links on the triangular 
lattice as described in Fig. 1(b). The model (3) is differ- 
ent from the known J1-J2-J3 model with all isotropic ex- 
change coupling previously studied in the literature 75-77 
because in our model the third term explicitly breaks the 
SU(2) spin symmetry which will have important influ- 
ences on the magnetic phases of the model, and because 
third-neighbor couplings are not considered. In the re- 
mainder of this paper we study the phase diagram of the 
Hamiltonian (3) by a combination of analytical methods 
and exact diagonalization on finite size clusters. 



second-neighbor spin-dependent hopping. The near- 
est neighbor hopping term gives a semi-metallic phase 
with Dirac nodes in the absence of second-neighbor hop- 
ping. However, the second-neighbor hopping term, t 1 — 
—fa + it" 'a 1 , is complex and spin dependent (originat- 
ing from the spin-orbit coupling), where 7 £ {x,y, z} 
indicates the red, green and blue second-neighbor links 
as shown in Fig. 1(b). a 1 and a are the usual 2x2 Pauli 
and identity matrices, respectively. The complex con- 
tribution of the second-neighbor hopping term (propor- 
tional to t") is a result of hopping via d-p-d ligands, while 
the real part is a result of direct overlap between d or- 
bitals, leading to the spin-independent amplitude f . 30 
Any nonzero value of t" immediately opens a gap in the 
spectrum and turns the semi-metallic phase to a Z2 topo- 
logical insulator phase. 30 

We include the Coulomb interactions by adding a Hub- 
bard term, 

n = H + UJ2n if n ii . (2) 

i 

The model (2) has been studied in the weak and inter- 
mediate interaction regime by use of slave-spin theory. 73 
The quantum spin Hall insulator found in the weak inter- 
action limit is unstable to a valence-bond solid (VBS), a 
very close relative of the expected AFM phase, for small 
spin-orbit coupling, and an exotic topological phase be- 
yond a critical spin-orbit coupling strength, t" . Since the 
slave-spin theory is a variational approach, it cannot cap- 
ture the possible magnetic phases which arise in the limit 
of strong Coulomb interaction. In order to address this 
weakness of the slave-spin theory, we analytically take 



III. SCHWINGER BOSON MEAN FIELD 
THEORY 

Recent numerical studies of the Hubbard model 
on the honeycomb lattice show a spin-gapped phase 
with no long-range correlations and no broken symme- 
tries at an intermediate regime of Coulomb interaction 
3.5<[//i<4.3. 38 ~ 42 In terms of a J x - J 2 Heisenberg model 
on honeycomb lattice, this spin liquid phase is located 
around J2/ Ji~0.06. 78 The existence of a spin disordered 
phase was further confirmed by other techniques, includ- 
ing functional renormalization group 76 , mean field and 
exact diagonalization, 77 but for higher values of J2I J\- 
One way wonder if the third term in Eq.(3), which is 
anisotropic and has some degree of frustration, may help 
stabilize the spin liquid phase. To answer this ques- 
tion, we employ Schwinger boson mean-field theory to 
investigate the stability of the spin liquid phase. We 
find the coupling J3 actually tends to stabilize a mag- 
netically ordered phase instead of disordered one. This 
technique has proven successful in incorporating quan- 
tum fluctuations 79 ' 80 and has beed applied to the frus- 
trated Heisenberg model on honeycomb, 81,82 kagome and 
triangular lattices. 83,84 

In the Schwinger boson approach the spin operators 
are replaced by two flavors of bosons at each site, 

Si = \h\Sh u (4) 

where = (6L,6jj,) are bosonic operators and a is the 
vector of Pauli matrices. For this to be a faithful repre- 
sentation at each site, the following constraint should be 
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imposed: 

fti=X>lX=2& (5) 

At the mean field level this constraint is imposed on av- 
erage, namely ifii) — 2S, which is taken into account by 
Lagrange multiplier A, taken to be independent of the 
site i. The exchange interactions can be written as fol- 
lows, which make the model suitable for constructing a 
mean-field theory, 

Si • Sj = X 0tij Xo,ij - A 0) yA ,y, 

Si • Sj — 2S?Sj = A Xti jA Xi ij — X x ,ijXx,ij, 

Si • Sj — 2SfSj — Ay^&yjj — X y ,ijXy,ij) 

S, • Sj - 2Sf 5? = A^.A^j - xl ijX z,ii, (6) 

where the A's and x's describe the paring and hopping 
of bosons. The full expressions are given in Appendix A. 

We use mean field-theory as a variational approach 
and decouple the above expressions in different channels. 
Hence, the bosonic Hamiltonian becomes 

H = Ji Xo,ijXo,ij ~ AS,ijAo,ij + h.c. 

<ij> 

Xo,ijXO,ij ^0,ij M ■ n * c - 

nnn 
nnn, a 

+\J2(^-2S)+E , (7) 

i 

where Eq is an energy constant. As usual, the minimiza- 
tion of the ground state energy with respect to mean-field 
parameters provides a set of equations which should be 
solved self-consistently. Although it is possible to solve 
these equations, in practice it is a formidable task to find 
the solution. We therefore use the pairings and hoppings 
as variational parameters which can be tuned by the cou- 
plings Ji, J2, J3. Moreover, we should note that since the 
representation in Eq.(4) has a U(l) gauge redundancy, 
the mean-field ansatz should be invariant under a com- 
bined physical symmetry group and gauge-group opera- 
tion which is called a projective symmetry group (PSG) 
operation. 85 In the PSG each physical symmetry is im- 
plemented followed by a particular gauge rotation such 
that the mean field ansatz is left invariant. We consider 
a uniform ansatz with zero-flux, 86 which would inspire 
a candidate for the short-range resonance valence bond 
(RVB) state. 87 The ansatz is defined as 



Xo,ij 


= Xi, 


A ,jj - 


A X) 


NN 




XO,ij 


= X2, 


A ,jj = 


A 2 , 


NNN 




Xcxjj 


= Xa, 




= A Q , 


NNN 


a — x, y, z 



Fourier transformed, the Hamiltonian can then be easily 
diagonalized as 

H=\^lH{k)*k + E», (9) 

k 




FIG. 2. (Color online) Schwinger Boson mean-field phase dia- 
gram of the model in Eq. (3) . The vertical axis is the average of 
boson density at each site, and the horizontal axis is the varia- 
tional parameter A2 /Ai which is a measure of the strength of 
second term in Hamiltonian. Different curves correspond to 
different values of J3 (A3 = 1). The dashed line indicates the 
physical value of spin, namely S = 1/2, which crosses the solid 
blue line in a very small window around A2/A1 = 0.5 indicat- 
ing the existence of a Z2 spin liquid phase. By increasing J3, 
the corresponding red (circle), square (green) and grey (star) 
curves evade a crossing of the dashed line making the spin 
liquid phase very unlikely. Note that the region below each 
curve is a Z2 spin liquid phase, and above it is a magnetically 
ordered phase of either commensurate or incommensurate. 



where = {b k AtbkAib k BthBibl kAt bl kM bl kBt bl kBi ), 
and 

(10) 

For full expressions of hu and A^ see Appendix B. 

With respect to the couplings Jj's the bosons may con- 
dense at some wavevector in the Brillouin zone. The 
zeros in the energy spectrum of the condensate is used 
to determine the magnetic ordering that develops in the 
system. 83 On the other hand, a gapped spectrum of 
bosons could signal the existence of a disordered phase. 
Hence, we can determine the phase boundary between 
condensed and uncondensed regions, i.e. ordered and 
disordered phases. For simplicity we take A a — A3. The 
phase diagram is shown in Fig. 2. Shown is a plot of 
the mean-value of the boson number at each site versus 
A2/A1 for different values of A3/A1. For A3 = there is 
a narrow region around A2/A1 « 0.5 86 which is crossed 
by the dashed line denoting the real value of S = 1/2. In 
this region the bosonic spectrum is gapped, thus the sys- 
tem is spin disordered. Because of the second neighbor 
pairing which "Higgs" the U(l) gauge bosons down to 
Z2, this spin liquid phase is stable 86 unlike the U(l) spin 
liquid phase found in Ref.[81]. Other parts of the phase 
diagram are magnetically ordered. At small values of J2 
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FIG. 3. (Color online) Top: Classical spin ordering of the 
Jx-Jz model (3) with J2=0. (a) Stripy order for J3 < x c i, 
(b) Neel order for x c i < J3 < x C 2, and (c) Spiral order for 
J3 > Xc2- Bottom: The phase diagram as a function of x = 
J3/J1 and the positions of phase transition points x c \ and 
Xci for the classical and quantum model. The results for the 
quantum model were obtained by exact diagonalization (ED) 
on finite size clusters. 



the ordered phase is a Neel phase where the bosons con- 
dense at the center of Brillouin zone, and for large values 
of J2 the condensations occur at finite momenta which 
would lead to an incommensurate magnetically ordered 
phase. 86 

Upon the inclusion of the third term J3, the spin liq- 
uid window disappears, and as seen in Fig. 2 the phase 
boundary is not crossed by the S = 1/2 line anymore. 
Therefore, it seems that the anisotropic J3-term strongly 
suppresses quantum fluctuations and favors a magneti- 
cally ordered state. This is in contrast to the J1-J2-J3 
model, where the third-neighbor term was shown to sta- 
bilize a spin liquid phase via the same Schwinger boson 
method. 82 (Recall that our J 3 is a second-neighbor cou- 
pling coming from t" in Eq.(l).) Despite the apparent 
breaking of SU(2) symmetry by the J3-term in Eq.(3), it 
has a "hidden" rotational invariance which could stabi- 
lize a magnetically ordered phase. We will elaborate on 
this point in the following section. 



IV. CANDIDATE MAGNETICALLY ORDERED 
PHASES 

The discussion in the preceding section showed that the 
J 3 term tends to stabilize ordered phases, and therefore 
disfavors spin liquids. In this section we discuss possible 
magnetically ordered phases at the limit of large spin val- 
ues, namely the classical orders. Ground-state configura- 
tions are readily obtained in some limiting cases. Since 
the model with J 3 = has been studied before, 88 here we 
focus on the J1-J3 model, and set J 2 — 0. We will discuss 
the effect of the J2 coupling in Sec. VII. In the limit of 
vanishing J 3 the magnetic phase is the usual Neel phase. 
On the other hand, in the limit of |J3|^>Ji, the model 
is decoupled to two trianglular latticies each governed by 
the J 3 term, namely H tri = - J 3 £) <ij><y S l • Sj - 2S]S 1 j , 



where, now < ij > stands for NN links on the triangu- 
lar lattice. This Hamiltonian, despite being obtained in 
an extreme limit, gives a proper low energy description 
of the cobaltates where the spin-orbit coupling is much 
stronger than the superexchange coupling. 89 The Hamil- 
tonian H tr i is not frustrated as can be seen by dividing 
the triangular lattice to four sublattices and performing 
the following local transformations on the four sublatti- 
cies shown in Fig. 1(b) by empty and blue circles, squares 
and triangles: 

Empty circle : S x = S x , S y = S y , S z = S z , 
Red triangle : S x = S x , S y = -S y , S z = -S z 7 
Green square : S x = -S x , S y = S y , S z = -S z , 
Blue circle : S x = -S x , S v = -S y , S z = S z . 

(11) 

In the transformed basis the Hamiltonian becomes 
fully isotropic, 



Htri — J3 Sj • Sj 



(12) 



Thus, the ground state will be the well known 120° order- 
ing on the triangular lattice for J 3 >0 and a fully ferro- 
magnetic state for Js<0. Transforming back to the orig- 
inal spins, we obtain stripy (or zig-zag) order for latter 
case and spiral order for the former case as shown, respec- 
tively, in Fig. 3(a) and Fig. 3(c). The stripy order was also 
observed in the Kitaev-Heisenberg model. 63 Note that 
the stripy order and spiral order result from the explicit 
SU(2) spin symmetry breaking. Therefore, any phase 
transition from the Neel phase to these phases should 
be related to a fully broken spin rotational symmetry, 
as such a phase transition is absent in the Kane-Mele- 
Hubbard model which preserves U(l) spin symmetry. 38 ' 90 
Having established the orderings in the extreme limits 
of the model, we can now estimate the classical transition 
points by comparing the ground state energy per spin: 

j? Eo 3 



E 



stripy 



E 



spiral 



2J 1 S 2 N ceU 
E 

~ 4J!S 2 N ce u 
' 2AJ!S 2 N cell 



3x, 



~2 X > 



(13) 

where x = J3/J1. For x < —0.25 the stripy phase is 
favored, for —0.25 < x < 2.75 the Neel phase, and for 
x > 2.75 the commensurate spiral phase is favored. In 
the next section we will show that the quantum effects 
shift the phase boundaries. 



V. EXACT DIAGONALIZATION ON FINITE 
SIZE LATTICES 

In this section we use Lanczos exact diagolalization 
(ED) on finite size lattices to locate the critical points 
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FIG. 4. (Color online) Upper panel: The ground state energy 
per site versus coupling J3 for different lattice sizes, from top 
to bottom: N=13 (grey), 17 (green), 19 (red) and 24 (blue). 
All clusters used open boundary conditions, and we set Ji=l 
and J2=0. Lower panel: The second derivative of the ground 
state energy versus coupling. Note the vertical axis is in a 
logarithmic scale. Two clear peaks become more pronounced 
as the system size increases. For the largest lattice size the 
left peak appears around J3=-0.62±0.01 and right one around 
J 3 =1.62±0.01. 



and identify different magnetic phases. We consider two- 
dimensional lattices with N=13, 17,19, 24 sites. For N=24 
we considered both periodic and open boundary condi- 
tions. Note that the periodic lattice with N=24 is the 
minimum lattice size whose extra frustration due to the 
boundary of the lattice is avoided. This cluster is shown 
in Fig. 1(b). The N=24 lattice avoids extra frustration 
from the periodicity of transformations in Eq.(ll), which 
is 2, because it is consistent with the periodicity of mag- 
netic ordering of Eq.(12) for J3 > 0, which is 3 for each 
direction on the triangular lattice. Naively, if we take a 
lattice based on the primitive unite vectors of the honey- 
comb lattice, it will be of size 72, which is not practical 
for ED. For sizes smaller than 24 the open boundary con- 
dition should be worked out. 

To find magnetic phases, set of ordered magnetic 
phases found classically in Fig. 3 serves as a useful guide. 
The idea is to define the corresponding order parame- 
ter for the stripy (zig-zag) and Neel phases. The un- 
frustrated model at J 3 =0 possesses antiferromagnetic 
long-range order characterized by a staggered magne- 
tization, with the staggered moment m=0.2677, which 
has been significantly reduced from classical moment by 
quantum fluctuations. 91 This has been verified by various 
methods including spin-wave theory, the coupled cluster 
method, exact diagonalization, series expansions around 
the Ising limit, tensor network studies, variational Monte 
Carlo and quantum Monte Carlo (QMC) simulations (see 
A.F.Albuquerque, et al 77 and references therien). For 
numerical purposes, a good order parameter of the Neel 
phase would be the staggered magnetization squared on 



entire lattice, 

" m+T) (? ( " 1)<Si ) • (14) 

where S = N/2. In the following, this staggered mag- 
netization is used to determine the stability of the Neel 
phase. At J\=0 the model is also unfrustrated thanks 
to the transformation in Eq.(ll). As pointed out in pre- 
ceding section for J 3 <0 (>0), the model reduces to fer- 
romagnetic (antiferromagntic) exchange coupling on two 
isolated triangular lattices, which can then described by 
the Hamiltonian in Eq.(12). Therefore, for J 3 <0 one can 
simply compile the following total ferromagnetic moment 
for rotated spins, 

m 2 stripy = s{s + 1) (E S ») ■ ( 15 ) 

These order parameters and their vanishing (to within fi- 
nite size limitations) describe the stability of the collmear 
ordered phases and their transition to a spiral phase. 

In Fig. 4, we plot ground sate energy (upper panel) 
and its second derivative versus coupling J 3 for differ- 
ent lattice sizes. The energy shows a smooth behavior 
with maximum around J3=-0.62±0.01 for the largest size 
N=24 with open boundary conditions. This maximum 
indicates that there is a phase transition at this coupling. 
The existence of a rather sharp peak appeared in the sec- 
ond derivative of the ground state energy (lower panel) 
clearly signals this phase transition. While it is rather 
hard to locate another phase transition just by looking 
at the behavior of energy, its second derivative shows the 
approximate location of the second-order transition to 
yet another phase. It occurs around J 3 =1.62±0.01 for 
N=24. One can also see that by increasing the size of 
the lattices the peaks in the second derivative of the en- 
ergy gets more pronounced. These critical values could 
be compared with the classical values in Fig. 3 which 
show that quantum fluctuations significantly shift them 
to lower values. 

The nature of magnetic phases around the critical 
points can be determined by examining the order param- 
eters in Eq.(14) and Eq.(15). Their behaviors are elab- 
orated for N=19 and N=24 in Fig. 5 and Fig. 6, respec- 
tively. For N=24 we presented the results for periodic 
boundary condition to avoid boundary effects. Besides 
the order parameters we also included the second deriva- 
tive of the energy to clearly show the singularity in the 
derivative of the energy coincides with the drop or onset 
of various order parameters. The nature of the ordered 
phases is already clear for the small size cluster N=19. 
For values of J 3 <-0.62 the dominant order is the stripy 
order (See Fig. 3(a)) which is steeply reduced across the 
critical point, where the Neel phase starts to develop over 
the intermediate range of coupling. The Neel order pa- 
rameter increases even beyond the isotropic J 3 =0 point 
owing to the ferromagnetic nature of the isotropic term in 
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FIG. 5. (Color online) Upper panel: Second derivative of 
the ground state energy per site versus coupling J3 for a lat- 
tice with N=19 and open boundary conditions. We set Ji=l 
and J2=0. Lower panel: Order parameters squared m 2 Neel 
(empty blue circles), m 2 stripy (green solid dots) and m 2 zig _ zag 
(squares). At the left critical point J3=-0.62 the stripy oder 
reduces and the Neel order starts to increase. Both orders 
decrease across the right critical point around J3 = 1.62 where 
a spiral order begins to develop. The vertical dashed lines are 
guides to the eyes. 
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FIG. 6. (Color online) The same quantities are plotted as in 
Fig. 5 but for N=24 with periodic boundary conditions. Note 
that the approximate positions of critical points have been 
changed. 



Eq.3. These two collinear phases then drop down across 
the second critical point, with the onset of spiral order. 

As far as the magnetic structure of the layered Iridate 
Na2lrC>3 is concerned, recent resonant x-ray scattering 92 
showed magnetic Brag peaks at wave vectors consistent 
with both stripy and zig-zag orders. However, first prin- 



ciples calculations 92 and very recent inelastic neutron 
scattering (INS) experiments 70,93 ' 94 showed that further- 
neighbor exchange interactions are strong, which in turn 
makes the zig-zag configuration the most likely magnetic 
order. 

Similar to the Heisenberg-Kitaev model, the model in 
Eq.(3) can not explain the zig-zag order as its ground 
state. However, we would like to argue that even in 
the context of the model in Eq.(3), we have some de- 
gree of zig-zag ordering. The argument traces back to 
the Hamiltonian of transformed spins on triangular lat- 
tice in Eq.(12) for J^<0 leading to ferromagnetic order 
of rotated spins. Transforming back to the original spins, 
the ordering on honeycomb lattice will be stripy or zig- 
zag depending on the sign of the NN coupling: stripy 
for Ji>0 and zig-zag for Ji<0. Often Ji> 0, so the 
former phase is energetically favored. But it seems quan- 
tum fluctuations make them close to each other in the 
strength of their respective order parameters. To elabo- 
rate this, in bottom panel of Fig. 5 we also plot the zig- 
zag order parameter of the ground state. It is clearly 
seen that, despite having higher energy than the stripy 
phase, the ground state posses a high degree of zigzag or- 
dering. Note that this holds at the level of an exchange 
model with up to second-neighbor coupling unlike the 
well known J\- J2-J3 75 ' 82 or Kitaev-Heisenberg- J2- J3 68 
model. The results for N=24 with periodic boundary 
conditions are shown in Fig. 6. The overall features of 
the plot is the same as Fig. 5 except the approximate po- 
sitions of the critical points have been displaced: The 
transition from the stripy to Neel and from Neel to the 
spiral phase would occur around J 3 =-0.4 and J3=1.3, re- 
spectively. Moreover, the value of zig-zag order is still 
high with a ratio m Z i g ^ zag / m str i P y— 0.62 at Jq=-1. 



VI. FIDELITY ANALYSIS OF THE PHASE 
TRANSITIONS 

In this section we use a quantum information theoretic 
tool called fidelity to analyze the quantum phase transi- 
tions described in the preceding sections. For pure states, 
say and \tp2), it simply measures the overlap between 
them as F = IV^)!, an d so is a measure of distin- 
guishability between the states. The states and j^) 
could be ground states of the Hamiltonian in Eq.(3) cor- 
responding to slightly different values of tuning param- 
eter, say J3 and J3+6, namely F( J 3 , S) = \ (ip( Jz)\4>{ Jz + 
6)) |, where S is a small quantity. Naturally for the same 
(orthogonal) states it will be unity (zero). While far 
away from the critical points the distinguishability is not 
significant, across a phase transition there is a drastic 
change in the fidelity as the ground states on different 
sides of critical points are totally different. Hence, the 
fidelity could be a strong signature of a quantum phase 
transition. 95 ' 96 Moreover, its scaling is intrinsically con- 
nected to derivatives of the ground state energy, and it 
was argued that the fidelity susceptibility (FS) defined 
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FIG. 7. (Color online) Fidelity (upper panel) and fidelity sus- 
ceptibility (bottom panel) versus J3 (with J2=0) for N=24 
with periodic boundary conditions and 5=0.01. The singular- 
ity in S(Ja) is much more pronounced than the singularity in 
the second derivative of ground state energy shown in Fig. 6. 
The dashed line indicates the location of critical point J3 = 1.3. 



as S(J 3 ) = a 5 2 F(J 3) 5)| 5=0 =21im^ is a more 

sensitive tool than the second derivative of the ground 
state energy to detect the critical points. 97 In particular, 
close to criticality, 

d 2 E JL_ 

dJl ~ A' 

S(Js) ~ ^, (16) 

where A is the energy gap. This relation explicitly in- 
dicates that a singularity in the second derivative of the 
energy will imply a singularity in the FS, and the sin- 
gularity in the FS is more pronounced close to phase 
transition. Moreover, as the definition of the fidelity or 
the FS does not rely on the notion of either a local order 
parameter or symmetry, it has been used to detect topo- 
logical phase transitions which evade a description based 
on a local order parameter. 98 ' 99 

In Fig. 7, we plot both fidelity (upper panel) and fi- 
delity susceptibility (bottom panel) for positive values 
of J3 and for a lattice size N=24 with periodic bound- 
ary conditions. Indeed, for Js<0 the phase transition 
is already clear from the derivative of energy as shown 
in Fig. 6. Of course, it is also clearly signaled in fidelity 
(not shown here). But the singularity in derivative of 
energy around J3=1.3 is rather weak and has a bump- 
like behavior which make the precise determination of 
the location of the critical point difficult. Fidelity, as ex- 
pected, exhibits a significant change in the ground state 
wavefunction across the critical point. Additionally, fi- 
delity susceptibility reveals a more pronounced singular- 
ity at the critical point J 3 =1.3 in contrast to the second 
derivative of energy in Fig. 6. 



VII. CONCLUSIONS AND FURTHER ASPECTS 

In this work we studied the strong interaction limit 
of the Hubbard model in Eq.(2) whose noninteracting 
ground state is a two-dimensional topological band in- 
sulator fully breaking the SU(2) spin symmetry, 30 which 
in turn makes the exchange couplings highly anisotropic 
and frustrated. Since the Schwinger boson mean-field 
theory shows that the anisotropic term favors magnetic 
ordered phases rather than a spin liquid phase which 
could arise in frustrated models, we first considered the 
J1-J3 model and set J2=0. The results of exact diagonal- 
ization studies on finite size lattices are summarized in 
Figs. 4, 5, 6, where it was shown that varying coupling J3 
(with J2=0) stabilizes the stripy, Neel, and spiral phases. 
For the largest lattice we considered, N=24, the criti- 
cal points are at J3—-OA and J3=1.3 corresponding to 
stripy-Neel and Neel-spiral transitions, respectively. We 
further exploited the fidelity susceptibility to locate the 
latter phase transition more precisely. The appearance 
of stripy and spiral orders, and the phase transitions out 
of a Neel phase are ascribed to the explicitly broken spin 
symmetry. 

These magnetic phase transitions can be used to shed 
light on correlation effects in two-dimensional topolog- 
ical band insulator. Indeed, the topological band in- 
sulator persists up to intermediate strengths of on site 
Hubbard interaction. However, the physics at intermedi- 
ate and strong interactions is significantly influenced by 
strength of the spin-orbit coupling, t" , which determines 
the gap of the non-interacting topological insulator. For 
intermediate values of the Hubbard interaction, small 
and large values of *" result in AFM (VBS) and QHS* 
phases, respectively. 73 While the former breaks the time 
reversal symmetry, the latter still preserves time reversal 
symmetry with protected collective edge modes and bulk 
topological degeneracy. Our work based on the exchange 
Hamiltonian in Eq.(3) further revealed that this topolog- 
ical phase eventually breaks down at strong interacting 
limit to a magnetically ordered phase with a spiral tex- 
ture. However, it appears that a phase transition as a 
function of spin-orbit coupling still persists in the large 
interaction limit. 

Thus far, we have considered the limit of J2—O. In the 
rest of this section we discuss the affect of this isotropic 
exchange coupling term and map out the full phase dia- 
gram for a lattice size N=19. The phase diagram in the 
J2-J3 plane is shown in Fig. 8. To obtain it, we used the 
knowledge of the phases appearing in various limits to 
distinguish different phases and the transitions between 
which are signaled by a singularity in the derivative of 
ground state energy. Although this lattice size is too 
small to locate the precise position of phase boundaries, 
we showed in preceding sections that enlarging the lattice 
just shifts the phase boundary somewhat. So, we believe 
that even systems of this small size can give a qualitative 
sense of the phase diagram of the model. 

On the J3=0 axis the existence of AFM, paramag- 
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Unlike the original Kitaev model 61 which is exactly solve- 
able due to trivalent structure of the honeycomb lattice 
giving rise to a bilinear representation of the Hamiltonian 
in term of Majorana fermions, an extension of the model 
to a triangular lattice is no longer trivial for spin- 1/2 de- 
grees of freedom. Even if we assume the "parallel" Kitaev 
models with isotropic couplings on the triangular lattice 
are spin disordered, the coupling between layers, though 
small, may stabilize a magnetically ordered phase. It is 
unclear what this order may be. We leave this question 
for future study. 
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FIG. 8. (Color online) Schematic phase diagram of the model 
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netic (PM) and spiral orders have been investigated in 
the literature. 76 ' 77 The only difference with the classical 
phase diagram 75 ' 88 is the appearance of a PM phase at 
intermediate values of J2- There are two spiral phases 
in phase diagram. The spiral phase on the Js=0 axis 
should be distinguished from the spiral phase arising at 
the large values of J3 (J2=0), which is why we showed it 
with a star. At the extreme limits where either J2^>J\ 
or J^>Ji the spiral* and spiral phases are adiabatically 
connected to the 120° ordering of, respectively, original 
and rotated spins (see Eq.(ll)) on decoupled triangular 
lattices. 

The phase diagram in Fig. 8 should be compared with 
the phase diagram reported for the isotropic J1-J2-J3 
model. 76 ' 77 Note that in the latter model J3 coupling 
stands for third-neighbor couplings on the honeycomb 
lattice. While this isotropic term stabilizes the AFM 
phase, the anisotropic coupling in Eq.(3) stabilizes the 
spiral order instead. However, in both models the PM 
phase persists up to intermediate coupling J3. There is 
still one phase labeled by a question mark in Fig. 8 which 
needs further study, being out of the scope of this work. 
However, we can argue which model may describe this 
undetermined phase. It appears this phase is stabilized 
when J2^Ji ft Ji- In this limit the second neighbor 
isotropic terms in Eq.(3) vanish and we will end up with 
a Kitaev-type model on two triangular lattices coupled 
antiferromagnetically through J\: 

///. = Ji Y, s i • s - + E s ? s l- 

<i,j> <Sij3>7 



Note added: Upon the completion of this work we 
became aware of the work in Ref. [90] which addresses 
the same model by use of the functional renormalization 
group. Part of our results, such as the identification of 
the Neel and spiral phases for J2=0, agrees except for the 
location of critical points. 

Appendix A: Expressions of Bosonic Pairing and 
Hopping Terms 

In this Appendix we provide some details of Schwinger- 
boson men-field theory. The exchange interaction in 
terms of pairing and hopping of bosons in Eq.(6) are de- 
scribed by the following expressions, 



XO.ij 


= 




A 0l ij 






Xx,ij 


= \(bl t b u 


+ b\ i b n ), 




= ^(Mit ^ 


- b H b il) 


Xy,ij 


= \(^b 3i 








biibji) 


Xz,ij 


= \{b\^ 


-blb 3l ), 




= \(b lt b u + 





(Al) 



Appendix B: Hopping and Pairing Matrices 

The hk and in Eq.(9) are hopping and pairing ma- 
trices, respectively, whose elements are, 
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hi 1 = J2X292 - >hXz cos(k ■ a 2 ) + A, /if = - J 3 [x x cos/c • (ai - a 2 ) + ix y sin(k ■ a{)}, h^ 3 = ~ JiXldU ^f = 0, 
h l 2 = J2X192 + 'hXz c °s(& • CL2) + A, hf = 0, /if = * Jixm, 

hi 3 = J2X292 - J3X.Z cos(k ■ a 2 ) + A, /if = -J 3 [x x cos A: • (ai - a 2 ) + ix y sm(k ■ ai)], 

/if = J2X292 + -hXz cos(/c • a 2 ) + A, (Bl) 

Af = Js[A x cos k ■ (ai — a 2 ) + A y cos(k ■ ai)], Af = iJ 2 A 2 {sm k ■ [ai — a 2 ) — sin(/c ■ ai) + sin(fc • 02)] + J 3 A Z cos(fc -02), 
Aj!, 4 = - JiAigi, A 21 = iJ 2 A 2 [— sinfc • (ai — a 2 ) + sin(fc • ai) — sin(fc • a 2 )] + JsA 2 cos(fc • a 2 ), 

Af = J 3 [-A,cosfc • (01 - a 2 ) + A y cos(fc • o x )], Af = -i^A^i, Af = -ijiA^J, 

Afe 3 = ^3[A X cos A: • (ai — a 2 ) + A y cos(/c ■ ai)], A^ 4 = iJ 2 A 2 [sin k ■ (ai — a 2 ) — sin(/c ■ ai) + sin(fc • a 2 )] + JsA x cos(fc -02), 
Afc 1 = ^JiAigi, A 4 . 3 = iJ 2 A 2 [- sin/c ■ (ai - a 2 ) + sin(A: • ai) - sin(fc • a 2 )] + J 3 A 2 cos(fc ■ a 2 ), 

Af = Jsl-A^cosfc • (ai - oa) + A„ cos(fc • 01)], A J, 3 = Af = Af = Af = 0, (B2) 
where 

,gi = 1 + e ~ ihai + e~ ik ' a2 , g 2 = cos(k ■ m) + cos(/c • a 2 ) + cos A: • (o a - a 2 ). (B3) 

I 
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